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Abstract 

We consider the interactions of two identical, orthogonally polarized vector solitons in a nonlinear 
optical fiber with two polarization directions, described by a coupled pair of nonlinear Schrodinger 
equations. We study a low-dimensional model system of Hamiltonian ordinary differential equations 
(ODEs) derived by Ueda and Kath and also studied by Tan and Yang. We derive a further simplified 
model which has similar dynamics but is more amenable to analysis. Sufficiently fast solitons move 
by each other without much interaction, but below a critical velocity the solitons may be captured. 
In certain bands of initial velocities the solitons are initially captured, but separate after passing 
each other twice, a phenomenon known as the two-bounce or two-pass resonance. We derive an 
analytic formula for the critical velocity. Using matched asymptotic expansions for separatrix 
crossing, we determine the location of these "resonance windows." Numerical simulations of the 
ODE models show they compare quite well with the asymptotic theory. 

PACS numbers: 42.65.Tg, 05.45.Yv, 42.81.Dp, 05.45a 



*Electronic address: goodman@njit.edu; URL: http://math.njit.edu/goodman/ 



^Electronic address: rhaberma@mail.smu.edu 



1 



I. INTRODUCTION 



Solitary waves are an important phenomenon in nonlinear physics and applied mathemat- 
ics. So itary waves have been studied in a diverse array of physical models including water 
wave, fl i& q _ e,ec tr on,c dev,ce S ( Jos epa S on JU n ctl o» s) Q aad co Smology fl Q 
One of the most important applications is to nonlinear optical communications where soli- 
tary waves have been proposed as information bits in optical fiber transmission systems []| 
and produced experimentally about 25 years ago [8||. Other solitary wave phenomena in 
nonlinear optics include gap solitons in Bragg gratings 0, and dispersion managed soli- 



tons 



which hold promise for eliminating the timing jitter associated with soliton 
transmission systems. 

A single solitary wave propagating through a uniform medium appears particle-like in its 
coherence and steady propagation. Of great interest are the interaction of multiple solitary 
waves and the behavior of solitary waves propagating through non-uniform media. Solitary 
waves of completely integrable equations are known as solitons, and their interactions can 
be described completely, using multiple-soliton formulas derived via the inverse scattering 
transform Q|. The infinite set of conservation laws in integrable systems severely con- 
strain the dynamics: collisions are elastic, and the solitons will re-emerge from a collision 
propagating with their initial amplitudes and speeds intact, although their positions will 
have undergone a finite jump. Solitary wave collisions in non-integrable wave equations, 
can usually not be found in closed form, and show a much richer variety of behaviors: the 
waves may attract or repel each other and, upon collision, the solitary waves may lose their 



coherence and break ap art. 



one another. 



14, 



16 



17 



merge 



18, 



19 



into a single localized structure, or even oscillate about 



20 



2l|. 



In a soliton-based communications system, the bits are represented by solitons. In the 
simplest scenario, the presence of a soliton in a given timing window codes a one, and 
its absence codes a zero. Collisions between solitons, coupled with random noise in fiber 
characteristics, can lead to large perturbations in the solitons polarizations and to timing 



jitter 



22| . A bit that arrives at the wrong time may be interpreted incorrectly by a receiver, 



as would a soliton that splits in half or two solitons that merge. Ueda and Kath show such 
behaviors are possible and cite several additional numerical studies of solitons collisions not 
included here. We an approach to the modeling and analysis of these phenomena that, while 



2 



highly idealized, leads to new insights into these collisions. 

Interacting pairs of solitary waves from several distinct (non-integrable) physical models 
have shown an interesting behavior in common. At high speeds, the solitary waves move 
right past each other, hardly interacting, while at speeds below some critical velocity, the 
solitary waves interact strongly and may merge into a single localized state. Interspersed 
among the initial velocities that lead to this capture are "resonance windows", for which 
the two waves approach each other, interact with each other for a finite time, and then 
move apart again; see the second and third graph in figure ^ This has been explored by 
Tan and Yang in a system of coupled nonlinear Schroding er ( CNLS) equations that model 
nonlinear propagation of light in birefringent optical fibers |l5l. Il6lll7|. and by Cambell and 
collaborators in kink-antikink interaction in the 4 equations and several other nonlinear 
Klein-Gordon models jisl . 19, 20, 21 1. These windows form a complicated fractal structure 



that has been described qualitatively and even quantitatively, but for which the underlying 
mechanism has been poorly understood. 

The same phenomenon was also observed by Fei, Kivshar, and Vazquez in the interac- 
tion of traveling kink solutions of the sine-Gordon and 4 equations with weak localized 



defects 



23, 



24, 



25]. Instead of two solitary waves merging, in this case the soliton could be 
captured, or pinned, at the location of the defect. Almost all of the described models have 
been studied using the so-called variational approximation, in which the complex dynamics 
of the full PDE are modeled by a small, finite-dimensional system of ordinary differential 
equations. 

The sine-Gordon equation with defect and the birefringent fiber-optic model discussed 
above feature a small parameter measuring the "non-integrability" of the system. In a re- 
cent publication, Goodman and Haberman 2(| , exploited this small parameter to construct 
approximate solutions to the system of ODEs for the sine-Gordon model derived in |24]. We 
calculated the critical velocity for defect-induced soliton capture via an energy calculation 
involving separatrix crossing, and the location of the resonance windows using a quantiza- 
tion condition that occurs in the asymptotic expansion. In the current paper, we apply the 
same method to derive similar quantitative features in Ueda and Kath's ODE model of soli- 
tary wave collision in coupled nonlinear Schrodinger equations, and to explain the structure 
underlying the fractal structure of resonance windows. 

In section[nJ we introduce the physical model — a coupled system of Nonlinear Schrodinger 
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equations — and describe previous results in which the "two pass resonance" phenomenon has 
been observed. In section IIII1 we introduce Ueda and Kath's finite-dimensional model sys- 
tem that captures the observed dynamics, and introduce a simplified model which partially 
linearizes the system and renders it amenable to our analysis. We show numerically that 
this simplification does not qualitatively alter the dynamics. In section I1V1 we set up the 
calculation as a singular perturbation problem, and describe the unperturbed dynamics. 
We determine the critical velocity by calculating the energy that is lost to vibrations as 
the solitons pass each other, employing a Melnikov integral. We generalize this calcula- 
tion slightly for subsequent interactions. In section El we construct approximate solutions 
using matched asymptotic approximations, incorporating the previously calculated energy 
changes. Section IV II contains a discussion of the differences between the original model and 
its simplification and presents a weakly nonlinear theory to account for them. We conclude 
in section IVIII with a physical summary and a discussion on the applicability of these results 
to other systems displaying similar behaviors. 



II. PHYSICAL PROBLEM AND PRIOR RESULTS 



Following the previously cited Q, |27|. we consider the model of polarized light propa- 
gation in a optical fiber, given by the system of coupled nonlinear Schrodinger equations 



id t A + dlA + (\A\ 2 + (3\B\ 2 ) A = 

(1) 

id t B + 2 B + (\B\ 2 + p\A\ 2 ) B = 0. 
This system replaces the more familiar scalar Schrodinger equation when polarization is 



taken into effect 



281 ] . The equations may be derived using the slowly varying envelope 
approximation to Maxwell's equations in an optical fiber waveguide. The variables A and 
B describe the envelopes of wave packets in the two polarization directions and f3 is the 
nonlinear cross-phase modulation (XPM) coefficient that arises due to cubic (x^) terms in 
the dielectric response of the glass. Here we use z as a space-like variable and t as a time-like 
variable. Of course, in the optics interpretation, the labels z and t are switched, as the signal 
is defined as a function of time at z = and the evolution occurs as the pulse moves down 
the length of the fiber. For mathematical simplicity, we will use t as the evolution variable. 
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Our interest is in the interaction of solitary waves in the above system. In the cases (3 = 
and (3 = 1, system (JTJ) is completely integrable (29, 30]. In the first case it reduces to a pair 
of uncoupled NLS equations; in the second it is known as the Manakov system. For other 
values of (3, the equations are not integrable. Of special interest is the case /3 = |, which 
corresponds to linear fiber birefringence. For very small values of (3, this system models light 
propagation in a two-mode optical fiber In the case (3 = 0, the equations are simply 
a pair of focusing nonlinear Schrodinger equations, with well known soliton solutions, first 
suggested as carriers of optical signals by Hasegawa and Tappert Q. When (3 takes any other 
value, the equations are non-integrable. Yang j^] studied these equations in great detail, 
enumerating several families of solitary waves and determining their stability. Of these, the 
only stable solitary waves come from a family of symmetric single-humped solutions. 

The simplest solutions of interest to (Q) consist of an exponentially-localized soliton in 
the first component, A, and zero in the second component, B, or vice versa. A single soliton 
propagates at constant speed with a fixed spatial profile. An important problem is the 
interaction of two such solitons upon collision, as interactions between two such solitons 
may lead to errors in a soliton-based transmission system. 

Tan and Yang numerically studied the interaction of two solitons initialized in orthogonal 
channels with identical amplitude, headed toward each other with exactly opposite initial 



velocities 



HQ 



la lia 113. For small values of (3 ~ 0.05, their simulations show that waves 



traveling above a critical velocity v c , the solitons pass by each other, losing a little bit of 
speed, but not otherwise showing a complicated interaction. At initial velocities below v c , 
the solitons capture each other and merge into a stationary state near their point of collision. 
See figure ^ which shows three graphs from 16J] in which the input velocity of the solitons 
is plotted on the x-axis, and the solutions followed until they separate, and their outgoing 
velocity plotted on the y-axis, and assigned the value zero if they merge. 

For somewhat larger values of (3 ~ 0.2, they find that in addition to the above behav- 
ior, that the capture region is interrupted by a sequence of "resonant reflection windows." 
Solitons with initial velocities in these resonance windows are reflected instead of being cap- 
tured. The numerical simulations show that the solitons pass each other once, undergo a 
finite number of width oscillations, then pass each other a second time. Thus they call this 
the "two-pass" resonance. 

For larger values of /3 ~ 0.6, they find not only reflection windows, but an intricate 
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fractal-like structure of both reflection and transmission windows. Certain portions of the 
structure, when properly scaled, look like copies, in come cases even reflected copies, of other 
portions of the structure, and such features are seen at many different scales. 

The two-bounce resonance in kink-antikink interactions was explained qualitatively in 



the first papers of the Campbell group [20, |2l|. As the kinks approach each other, they 



begin to interact, and, at time ti, and energy is transferred into a secondary mode of 
vibration, with some characteristic frequency to. If the initial velocity is below a critical 
value, the kinks no longer have enough energy to escape each other's orbit, and turn around 
to interact a second time ti- They show numerically that a resonant reflection occurs if 
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ti — t\ ~ lunjuj + 5. The parameter 5 is found by a least squares fit with numerical 
data. This relation is used to estimate the resonant initial velocities. This reasoning has 



subsequently been adapted in studies of sine-Gordon kink-defect interactions 



vector soliton collisions 



15 



16 



24, 



251 ] and of 



which are the focus of this paper. 



III. THE MODEL EQUATIONS 

In order to gain further insight into the resonance phenomenon, Tan and Yang examine a 



model system derived by Ueda and Kath [2j| using the variational method. In the variational 
method, the solution is assumed to take a certain functional form A(p(t)), B(p(t)), dependent 
on parameters p(t) that are allowed to vary as a function of time. This ansatz is then 
substituted into the Lagrangian functional for the PDE, which is integrated in space to 
yield a finite-dimensional effective Lagrangian, 



/oo 
C(A,A*,B,B*)dz, 
-oo 



whose Euler-Lagrange equations describe the evolution of the time-dependent parame- 
ters. Equation (JTJ) has Lagrangian density. 

C = i{AA* z - A Z A*) + i{BB* z - B Z B*) + (\A t \ 2 - |A| 4 ) + ( | ^ £ | 2 - |5| 4 ) - 2(3\A\ 2 \B\ 2 (2) 

Many examples using this method for PDE's arising as Euler-Lagrange equations are given 
in a recent review by Malomed 1331 . 

n 

Following [2j|, we take an ansatz corresponding to two solitons at distance 2X of equal 
magnitude heading toward each other with equal speed, 

z — X { b \ 

A = 7] sech exp i I viz — X) H (z — X) 2 + a ) , 

w \ 2w J 

z + X ( b x 

B = -nsech expi I — viz + X) H (z + X) 2 + a 

w \ 2w 

where 77, X, w, v, b, and a are time-dependent parameters for the amplitude, position, 
width, velocity, chirp, and phase, whose evolution remains to be determined. The variational 
procedure yields a conserved quantity K = r] 2 w, related to the conservation of the L 2 norm 
in CNLS, as well as the relations ^ = v and % = b. The evolution is described by the 
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Euler-Lagrange equations 

d 2 X 16K0 d 



F(a) (4a) 



dt 2 w 2 da 

10 n K-3(3K^-(aF(a))} (11,) 



dt 2 -k 2 w 2 \w da 
where a = X/w and the potential and coupling terms are given by 



x cosh x — sinh x 

F (x) = —3 • (5) 

sinh x 



Note that F, actually —F, is a potential term, not a force. We keep this notation for 
continuity with previous studies. 

Numerical simulations show that for small /3, a solution to with an initial condition 
of the form given in ansatz © will remain close to that form, i. e. the solution will continue 
to consist of two nearly orthogonally polarized solitons, at least until they merge into a 
single bound state. Using the symmetries of equation (jlj), we may set K = 1 without loss of 
generality. Equivalently, the PDE symmetry may be used to set K — 1 in the ansatz used 
by the variational method. 

These equations display the two-bounce resonance phenomenon, as shown by Tan and 
Yang. Consider the initial value problem, with "initial" conditions describing the behavior 
as t — > — <x>. 

dX dw 
X -> -oo; — — >• v in > 0; w — ► 1; — -> 
at dt 

This does not strictly determine a unique solution, since the solution is invariant to time 
translation. We plot t> out as a function of v- m with (3 = 0.05 in figure CD These and all other 
ODE simulations were performed using routines from ODEPACK |34|. Compare this figure 
to the three plots of figure ^ There are key similarities and differences between this graph 
and the exit velocity graphs of the full PDE simulations. The critical velocity in this figure is 
about v c = 0.19, close to the value v c = 0.1715 found in Q]. A noteworthy difference is the 
complex behavior of solutions with initial velocity below v c — no such behavior, not even the 
two-pass windows, was seen in the very careful simulations of Tan and Yang. This should 
not be surprising, as system (jlj) is Hamiltonian, and the set of initial conditions leading 
to unbounded trajectories in backwards time and bounded trajectories in forward time has 
measure zero, by reasoning similar to Poincare recurrence, as shown in Proposition 1 of [35]. 
Localized solutions to (0) may lose energy to radiation modes, a dissipation mechanism not 
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FIG. 2: The input vs. output velocity of a pair of orthogonally polarized solitons with j3 = 0.05 



present in the ODE model. As a further result of the dissipation, the output speeds of the 
reflected solutions are much smaller than the input speeds in the PDE solutions, whereas at 
the very center of the ODE windows, the output speed exactly matches the input speed. A 
more interesting difference can be seen in the presence of the wide reflection windows, which 
were not found in the PDE simulations with this value of (3, summarized in figure ^ 

In figure H3 the exit velocity graph of (@J) shows that even at /3 — 0.2, the ODE dynamics 
display a complex fractal-like structure in addition to the reflection windows, which are not 
seen in the PDE dynamics for such small values of f3. The numerical value of the crititcal 
velocity is v c = 0.86, close to the value v c = 0.936 found in [l6^ . 

The numerical solutions of (jlj) qualitatively explain the resonance windows. In figure 
the w(t) components of the solutions with initial velocity v at the center of the first two 
resonance windows (actually the points tangents to the line t> ou t = — v m)- In the leftmost 
window, the oscillator w(t) is excited, oscillates about 5 times and then is de-excited. In 
the next window, w(t) oscillates 6 times. In each of the successive windows, w(t) oscillates 
one more time before it is extinguished. We will refer to the first window as the 2-5 window 
and the second window as the 2-6 resonance window. Recall that no such windows have 
been found in the PDE dynamics for this value of /3, but such windows have been found the 
in the ODE dynamics for all values of j3. Tan and Yang demonstrated a width oscillation 
in the PDE solutions in analogy with that shown here. The minimum value of n in the 
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FIG. 3: Left: The exit velocity graph for equation (j3J) for (5 = 0.2, showing reflection windows and 
a variety of more complex fractal-like structures. Right: The same figure with all but the main 
resonant reflections removed. 
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FIG. 4: Plots of the w{t) component of (@J) with initial velocity v = 0.09988 (top) and v = 0.13464 
(bottom) and (3 = 0.05, showing the 2-5 and 2-6 resonances. 

2-n resonance decreases with increasing j3. There does exist a 2-1 resonance with velocity 
v = 0.649 in the ODE dynamics shown in figure H3 while the first resonance window found 
in the PDE simulations is the 2-2 resonance at about v = 0.9 in figure ^ 
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A. A further simplified model 



The model (J3J bears a striking resemblance to the system derived in 24| t° study the 
two-pass resonance in the sine-Gordon equation with defect, and analyzed in 2f|. In that 
case, however, the situation is much simpler: the term equivalent to w in (@J) occurs only 
linearly, and the potential and coupling terms, equivalent to F(X/w) here, are functions 
of X alone. This allows us to solve the analog of (j4b|) by variation of parameters to solve 
for this term and then insert it into the equivalent of equation (|4aJ). a critical step in our 
analysis. In our numerically-computed solutions displaying the two-bounce resonance for 
small values of (3, the width w undergoes only a small oscillation about its initial width 
w — 1. Therefore, we may partially linearize system (jl]), which allows us to proceed in the 
same manner as we have for the sine-Gordon system. We find reasonable agreement, with a 
few notable differences, between the two ODE systems. We will discuss the linearized theory 
first and then discuss corrections due to the nonlinearity. 

Allowing w = 1 + W, where W is considered small, expand all the terms in W, and keep 
only leading-order terms. We arrive at the reduced system: 

^i = 16(3(F'(X)+G'(X)W); (6a) 
d 2 W 16 TI7 48/3 

-^ + -^=JG(X), (6b) 

where 

G(X) = -(XF(X))'. 

Figure |3] shows that this simplified equation gives an accurate estimate of the critical velocity 
for small values of (3 based on numerical simulation. Figure |H1 shows the equivalent of figure El 
with the same value of f3 = 0.05 for the simplified equations. It shows that the qualitative 
picture, chaotic scattering interupted by resonance windows for v < v c , is the same, while the 
actual location of those windows varies greatly. In the case (3 = 0.05, the simplified equation 
has a 2-4 window, while the full equation's first resonance is 2-5. As f3 was decreased further, 
the agreement between the two systems improved. We rescale time by allowing t — > ^yffit, 
transforming equations (@J) to 

X = F\X) + G'(X)W; (7a) 
W + X 2 W = ^G(X), (7b) 
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FIG. 5: The critical velocity for capture as a function of the coupling (3 for the fully nonlinear 
system (J3J) (solid) and the simplified system (jSJ) (dashed). 
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FIG. 6: The exit velocity graph for the simplified system (jEJ), showing qualitative agreement with 
figure [21 



with fast frequency A given by 

A = (8) 

The dot notation will be used for derivatives with respect to the scaled time. The conditions 
in backward time as t — > — oo become: 

X -oo; X -> V m > 0; W -> W -> 0. (9) 
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We will use capital V to represent velocities in the scaled time t and lower-case v for velocities 
in the physical time. 



IV. DETERMINATION OF ENERGY CHANGE AND CRITICAL VELOCITY 
A. Setup of Melnikov Integral for AE 

First, note that if W is held equal to zero, equation (|7aj) has the phase space shown in 
figure showing three distinct types of orbits: closed orbits, corresponding to a pair of soli- 
tons bound together as a breather, unbounded orbits, corresponding to two solitons passing 
each other by, and orbits heteroclinic to degenerate saddle points at (X, X) = (±00, 0)— 
separatrices — that form a boundary between the two regimes. These orbits correspond to 
level sets, where the energy 

E = l -X 2 - F(X) (10) 

is negative, positive, or zero, respectively. As W is allowed to vary, solutions may cross the 
separatrices. We will show below that W remains 0(y/j3) by variation of parameters ()13|) 
below and, thus, that perturbation methods are applicable. 

We wish to asymptotically analyze orbits near the separatrix (see figure |7J), since two 
solitons are initially captured when they cross the separatrix and are reflected or transmitted 
when they cross it a second time and may escape. We first determine the energy loss as a 
soliton goes from X = —00 to X = +00 by computing an energy integral called a Melnikov 
integral [36] . A Melnikov integral is a perturbative device for measuring the change of energy 
in a given system. It is simply the integral of the time rate-of-change of the energy along 
some trajectory in the unperturbed problem. A zero of the Melnikov integral is commonly 
a necessary condition for chaos in low- dimensional dynamical systems j^. In our case, we 
simply wish to calculate a change in energy. 

The calculation has been simplified significantly from that given in [2^ , in a manner that 
yields additional insight into the form of the energy loss. In particular, we do not need to 
keep track of whether certain functions possess even or odd symmetry, and we find in an 
elementary way that the change of energy is negative. First, we note that the separatrix is 
given by the level set E = 0, therefore, along the separatrix, equation (fTTH) may be solved 
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FIG. 7: The X phase plane, showing trapped (dashed), untrapped (dash-dot), and separatrix (thin 
solid) orbits, corresponding to level sets of lfTU|). Superimposed is the X — X projection of the 2-6 
resonant solution to the fully nonlinear equations @ with (3 = 0.05 (thick solid line). 

for X(t), giving 

^ = vmxV)- (ii) 

Given the function F in ((HJ), it is not possible to find the separatrix orbit Xs(t) in closed 
form. The time-dependent energy exactly satisfies the differential equation 

^ = (X - F(X)) X = XG'(X)W = (^G(X(t))J W, (12) 

where we have used equations (|7aj). We approximate the change in energy for one nearly het- 
eroclinic orbit along the separatrix (from one saddle at infinity to the next saddle approach) 
by approximating X(t) in (|12|) with the known separatrix solution Xs(t). We integrate f)12|) 
along the length of the orbit and integrate by parts to find the total change in energy: 

* E= m G{xs{t)) ) wdt 

where we have integrated by parts. Given the initial condition Q, with V — for the 
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separatrix case, we may solve equation (|6b|) for W(t) using variation of parameters: 



Wit) = -^cosAt [ G(X s (r)) sin Xrdr + -^-sinXt f G(X s (t)) cos At dr, (13) 

7T 2 A ./_ ^ 7C 2 X 



(again approximating X(t) by X s (t)) and 
dW[t) 



3 /"* 3 f* 

— sin At / G(X s (r)) sin Ar dr H — -cos At / G(X s (r)) cos Ar dr 

M 7T 2 TV 2 

Setting I s (t) = J_ G(Xs(t)) sin Xrdr and I c (t) = G(Xg(r)) cos Xrdr, we find that 
AE = -- I s (t)I' s (t)dt - - / / c (*)/#)dt 



3 r(/ s 2 (oo) + / c 2 (oo)). 



(14) 



2tt 2 

This may be integrated by a standard substitution to yield 

AE = "i((/ G ( X s(r))sinAr dr^j + ^ G(X s (r)) cos Ar dr 
G(X s (r))e iAr dr 



3 



2tt 2 



(15) 



Thus, the problem is reduced to to calculating the integral ()15j) . In fact, because in this case 
G(Xs(t)) is an even function, 

AE = -^/c(oo) 2 . (16) 

Note that this shows the change in energy is generically negative when we assume W — > 
as t — > — oo. In fact, it must be negative, as the system conserves an energy that is 
positive-definite as \X\ — > oo, and no energy resides in the width oscillation initially. Using 
AE = -f , we find 

_V3 T 

7T 

where J Cj00 = J c (oo). The integral in equation (fTT)|) may be solved numerically by converting 
it into a differential equation, which may be integrated simultaneously with equation (|lljl. 
Alternatively, we derive closed-form approximations to AE v c in section flV CI below using 
complex analysis. 
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FIG. 8: The critical velocity of the ODE system © computed via direct numerical simulation 
(solid line), the first-order asymptotic approximation (dots) and the second-order asymptotic ap- 
proximation (dashed line) from (j24j) . and via numerical evaluation of integral (|15[) . 

B. A Generalizaton of the calculation 

Next, we briefly mention two generalizations of the Melnikov calculation above that will 
be useful later. First, suppose that instead of approaching zero as t — > — oo, 

W~ 4vWVsinA(f-0), 

71 A 

Then the change of energy will be given by 

A£ = ^(-i + WcosA<^. 

As X traverses the heteroclinic orbit in the reverse direction, the sign of G'(Xs(t)) in ()12|) 
is reversed, which leads to: 

AE = ^(~-WoobxA. (17) 

For a resonance to occur, the change of energy calculated in the first Melnikov integral must 
cancel with the energy jump on the return trip. Assume the forward heteroclinic orbit has 
"symmetry time" ti at which X = 0, with symmetry time t 2 on the return trip. Then, 
by (0, as t — > oo on the forward heteroclinic orbit, 

w(t) ~ 4r J ^ sinA ( t - t i)- ( 18 ) 
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For an exact resonance to occur, the energy change along the two heteroclinics must cancel, 
leading to the condition 



+ AE 2 = -^SS + (-1 - C os X(t 2 - tx) ) = 0, 

2n z tt z \ 2 



obtained by combining (|15|) and (|17|) with W = 1 and <fi = t 2 — t%- Thus cos A(t2 — ^i) = — 1, 



(2n+l)vr . , 

t 2 -t l = K - j-L-. (19) 



or 



This differs from the equivalent resonance condition in |26(, in which t 2 — ti = 2Tin/\. The 
difference arises because in that system the equivalent term to G(X) was an odd function, 
whereas here G is even. 

Many analyses of two-bounce and two-pass resonance phenomena have been based on the 
assumption that t 2 —t\— + 5 for some undetermined 6, a phase shift that accounts for 
unidentified physical processes that have not been modeled. Equation ()19|) shows that in 
this case 5 = f . It is worth computing a linear fit of t 2 — 1\ vs. n for comparison with earlier 
studies and, we will see, for comparison with the analogous result for the fully nonlinear 
ODE (@J). At p = 0.05, we find the linear fit t 2 -h= 4.935 (n + |) - 0.011 and for (3 = 0.2, 
we find t 2 -h w 4.931 (n + |) + 0.139, whereas ^ ps 4.9348. Therefore, we see that to 
leading order, relation (fjTij) holds, and that the agreement improves with decreasing j3. 

C. Evaluation of critical velocity using complex analysis 

Since A is large and G(X) is analytic, the calculated change of energy (|15|) is exponentially 
small. In a calculation for a similar system, we were able to calculate the analogous integrals 
explicitly because Xs(t) was known in a simple closed form J^. In the present case, we 
are forced instead to expand the integrand of (JT5j) about a certain (branch) pole. Given 
the form of the potential F has a pole whenever sinhX = and the numerator of F is 
nonzero or has a zero of order less than 3. The nearest pole to X = occurs at X = in. 
Along the separatrix 

— = y/2F{X), (20) 
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so we let t* be chosen such that X(t*) = given the initial condition X(0) = 0. This gives 
the formula: 



dt = 
'0 io 

We expand (J2*U|) about X = ?7r and t = t* and find 



/ / / r// / J—L—dX=-^=f \ ciy« 2.10392 z. (21] 

' ' W 2F(X) v^y Vsiny-i/cosy 1 y 



z — (X - in fl 2 + 0{{X-iufl 2 )^ dX= dt 

( Z 1)^V2 (X _ i7r)5/2 + 0((x _ m) n/ 2) = t _ f 

5V7T 

which may be inverted to form 

X — Z7T — (-l) 1 /10 7r V5 2 - 1 /5 5 -2/B(t _ f )2/5 + 0((t _ f )8/5) _ 

Based on the expansion 

F(X) = iir(X - in)- 3 + 0(1) 
we compute the two leading order terms of the integrand of ([To]) 

/_i\3/5 6/5 2 4/5q / i\l/5 2/5o8/5 

G(x s (t)) = ( 1} J /5 d (t - *r 8/5 + ( j J /5 (t - *r 6/5 + oat - tr 2/ % 

Therefore (fTo^l involves integrals of the type 

/oo 
e iXt {t-iTfdt 
-oo 

with A > 0, T > 0, and p < 0. Here z'T is the branch pole, from which a branch line extends 
vertically to zoo. By a shift of contour and a change of variables to z — i\(t — iT), this can 
be replaced by an integral over the Hankel contour 7, which starts at —00 below the real 
axis, circles zero once in the positive direction, and returns to —00 along (and above) the 
real axis 38^ . 

I(\,T,p) = { -^-e~ XT j e^dz 

which forms part of a familiar representation of Euler's gamma function and yields the 
exponentially small term 

I(X,T,p) = (-l)- p/2 2sin((p+ l)7r)r(p+l)e~ AT A- (p+1) . (22) 
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Using (|22j) and standard trigonometric and gamma function identities, we evaluate the 
integral in (JT5j) 

f * - (^T^r (?) sm + g) sin ^ + O(A-) 

(23) 

As the integrand is real, we choose the branch (— 1) 1//5 = — 1 above. Using that AE = y 
and the scaling relation given before (J7J), we arrive at the expansion for the critical velocity 
in physical variables: 

Vc = ^l e ~T/^ ^( 2/ 5) a/S i/5 _ e(A/5)a 2 5 + ...), (24) 

where 9(x) = sin7rxr(x) and a = 7r~ 2 / 5 2 4//5 5 2 / 5 . using equations (jSJ), fjiup . as well as the 
two integrals above. Figure |H1 shows that the critical velocity is poorly predicted by the first 
term in this series, but well-predicted up to about (3 — 0.1 when the second term is added. 
The series expansion of the integrand of (JT5|l about t = t* contains one more integrable 
term which does not lead to a visible improvement of the approximation to v c . In order to 
improve the approximation, one would have to calculate expansions about the additional 
singularities of G(Xg{t)) further off the imaginary t-axis. 

V. MATCHED ASYMPTOTIC CONSTRUCTION OF SOLUTIONS 
A. The expansion framework 

If Un > Vc, then X remains positive for all time and X — > +oo monotonically. We can 
call this a one-pass transmitted solution. A "pass" will occur each time X = 0, when the 
two solitons pass each other and energy is transferred between the translation and vibration 
modes. If V- m < V c , then the energy is negative after one pass, and the solitons reverse 
direction, setting up the second pass. On the first pass, the change of energy was shown 
in (JT5j) to be negative, but on subsequent passes, it may take either sign, by (JT7J). On the 
second, and subsequent, passes the solitons may escape if the energy is positive, or may be 
reversed again. We will focus primarily on the case that the solitons interact twice before 
escaping. 

Following j3| , we construct 2-pass solutions by a matched asymptotic expansion. The so- 
lution consists of sequences of nearly heteroclinic orbits connected to near saddle approaches 
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at X = ±00. The change in energy from one saddle approach to the next is approximated 
by the Melnikov integral calculated in section ITVl The 2-bounce solution can be constructed 
from the following 5 pieces: 

1. A near saddle approach to X = — 00 with energy E = such that X -> V in < V c , 
as £ \ —00 ; 

2. a heteroclinic orbit with dX/dt > such that X(t\) = 0, with energy change AE\ 
given by (fTtjj); 

3. a near saddle approach to X = +00 with negative energy E = — |M 2 , such that X 
achieves its maximum at £ = £*; 

4. a heteroclinic orbit with dX/dt < such that AT (£2) = 0, with energy change AE2 
given by fTTj) ; 

5. and a near saddle approach to X = —00 with positive energy E = \V^ nt , such X — > 
-Kut, as t Z 1 00. 

The times £1, £2, and £*, as well as the energy levels, remain to be determined below. In the 
language of matched asymptotics, the approximations at steps 1, 3, and 5 are the "outer 
solutions" and steps 2 and 4 are the "inner solutions." 

A comment about the last step is in order. For general initial velocity V- ini the energy 
at step 5 will not match the energy at step 1. If these two energies match exactly, then 
we say the solution is a 2-pass resonance. If the energy at step 5 is positive but less than 
VJn/2, then the solution is in the 2-pass window, and may be called an incomplete resonance. 
Physically, the solitons reflect off each other, but with reduced speed and with significant 
energy remaining in their width oscillation. The outer edges of the window will be given 
by velocities where the energy at step 5 is identically zero. This defines the width of the 
windows. If at step 5, the energy is instead negative, then the solution remains trapped 
for another step, alternating between negative energy near-saddle approaches to X = —00 
and X = 00 until enough energy is returned to X such that E = V^ ut /2, and X — * ±00. 
Non-resonant solutions and higher resonances are explained in section IV El 

For the simpler sine-Gordon system, we wrote down a general asymptotic formula for 
n-pass solutions, calculated the location of 3-pass windows, and calculated the widths of 
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the 2-bounce windows [26]. Analogous results are possible in the present situation and are 



discussed below in section fVEl although in less detail than in the previous paper. 



,40]. 



26, 39, 401. The heteroclinic 



We use the method of matched asymptotic expansions, as in 
orbits along the separatrix are matched (forward and backward in time) to the finite time 
singularities associated with the near-saddle approaches. We will not make use of the two 
positive energy expansions, so we will not compute them. They enter the analysis when the 
energy change calculated over the heteroclinic orbits in the above section will then be used 
to connect the positive and negative energy expansions. 



B. Asymptotic description of heteroclinic orbit for large X 

We first construct an expansion of the "inner solution," given by the heteroclinic orbit. 
Along the heteroclinic orbit, X 2 /2 = F(X). Setting X = at t = t\, the trajectory is given 
as the solution to 

f l dt' - I* dX ' [ X ° dX ' I* dX ' (25) 

Jti Jo v /2F(x 7 ) Jo V^vn Jx y f mx r ) 

for an arbitrary 0(1) constant X < X. The first integral is 0(1) and will be asymptotically 
dominated by the second. For X > 1, we may approximate the potential by 

F(X) ~ 4(X - l)e- 2X . (26) 

If we let Z = X - 1, then 

F(X) ~ (4e- 2 )Ze- 2Z , (27) 
which we substitute into (|25|) and get 

e f e z 



t-t 1 = — / —dZ + 0(l). 

VsJ Vz 

We make the substitution Z = C + y, where C> 1 and y = 0(1), as motivated in the next 
section, and expand the integral in powers of O -1 , yielding 



e 1+c ( „, 1 „ / 1 

/8C V~ ' 20 v ~ ' " \C 2 



t -^ = ^\ ey + ^-y)e y + 0[ — )) +0(1) (28) 
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C. Asymptotic description of the saddle approach near X = oo 



As the width perturbation W(t) remains small, we may construct approximate solutions 
from solutions to equation (|7ajl with W — 0. Since F(X) — » as |X| — > oo, (oo, 0) is a 
degenerate fixed point and is of saddle-type. First we compute the near-saddle approach at 
X = +oo, under the approximation that the solution has small constant energy given by 

2 

with M <§C 1. So that, using (J27J) . the near-saddle expansion satisfies: 

1 . M 2 
-Z 2 -Ae~ 2 Ze- 2Z = - — . 

2 2 

We make the expansion Z = C + y, where C ^> 1 is determined in (|29jl and show y = 0(1). 
Then 

^ 2 -4e- 2 (0 + 2/ )e- 2C e-^ = -^. 

We define C by 

80e" 2 e- 2C7 = M 2 (29) 
and let T = Mt, which gives the simplified equation 

_ u m r\ P -*v = _i. ( 3o) 



c/ry \ c 

D. The matching procedure for near and exact 2-pass resonances 

The simplest asymptotic approximation 

Ignoring the term ^ in (|30p. as C ^> 1, the near saddle approach takes the form 

e y = cos (T — T*) 

where T* is the "center time" at which the near saddle approach comes closest to the de- 
generate saddle at (X,X) = (oo, 0). This has finite-time singularities forward and backward 
in time. Backwards in time, this is singular as T — T* \ — |, and may be asymptotically 
expanded as 

e v = T - T* + |. (31) 
22 



For large X, the heteroclinic orbit in (|28|) is asymptotically approximated using (J2HJ), to 
leading order in C^ 1 , by 

e v =T - Ti 

with a similar expression for T 2 — T along the return heteroclinic. The algebraic growth of 
the heteroclinic orbit matches to the finite time singularities of the near-saddle approach 
only if 

7T 

T* - Ti = -, 
2' 

which, combined with a similar relation for T 2 — T*, yields T 2 — 7\ = tt as in our analysis of 
the sine-Gordon model, or, in the unsealed time variable 

which shows the energy dependence of the period. The energies at steps 1 and 3 are related 
by V -f + AE X = 

Exact Resonance Condition: 

To this point, the calculation has been valid for general Vj n < V c . We now specialize to 
the case of exact 2-pass resonance. In a resonant solution the second energy jump must 
balance the first AE 2 = —AEi, a condition for which is given in ([19)1 . implying X/(2n + 1). 
Thus the resonant initial velocity V n solves 

V 2 M 2 

— + AE 1 = (33) 

2 2 K ' 



Solving this for V n , using that = —V 2 /2, 



(34) 



7T 2 (2n+ l) 2 ' 

where scaling (JSJ) has been used to convert this result back to the physical variables. Also 
nmm{/3) is given by the smallest integer n that makes v n a real number in equation (|34j) . In 
figure El we see that this does a relatively poor job at predicting resonant velocities. 

An Improved Approximation for Resonance 



The above calculation showed clearly the procedure used to find the resonant velocities. 
Here it is refined slightly to improve its accuracy. After some rearrangement, equation (j3*Uj) 
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becomes 

dV = dt (35) 

(l + _g_) e -* - 1 

Scaling the time variable as before, we expand this integral in powers of J- and keep the 
first two terms. This gives 

dy 1 ye~ 2y dy 

T^TT + 2C n (-1 + e^yfl* 

for the near saddle approach, which has solution: 



dT (36) 



v ■ ^ T-T* \ 
e y — sin — — ; = cos j— . 37 

As T — T* \ ^1 + f on the right-hand side, we find y — > — oo. Rearranging the 
expansion of the heteroclinic orbit, we find (J28|) 

1 \ y T-T\ 



2C n 



Matching these two approximations gives T2 — T\ = ^1 + 7^-^ tt or 



1+ 2^ K 

a more accurate approximation to the period. Combining this with (|19|) yields: 



1 + 2k) 7T (2n + l)n 



M n A 

This is still not a closed equation as C n has yet to be specified. We may eliminate M n from 
this equation and ()29j1 to obtain an implicit relation that defines C n : 

32vr 2 /3(2n + l) 2 Cl - (2C n + \f e 2Cn+2 = 0. (39) 

This has exactly two positive roots as long as 27r 2 e _2 /3(2n + l) 2 > ^|^ e , with the larger 
root relevant. Thus, we come to the revised estimate of the resonant velocities 



2 16 I 1 



(40) 



Vn ~\ V " vr 2 (2n + l) 2 ' 
Figure El shows that this does better than our first estimate. In a similar computation, we 
found that this analysis in a neighborhood of X = oo was enough to determine the resonant 
velocities. We find in the next section that we can do better with a numerical criterion based 
on (HU). 
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A Numerical Condition for Resonance 



In all situations where heteroclinic or homoclinic orbits are matched to near-saddle ap- 
proaches t 2 — ti equals half the period of (J7a|) with W set to zero, 



X - F(X) = 0. (41) 

Let P(E) be the period of the closed orbit of ()41|1 with energy E < 0. We may solve for 
P{E) by evaluating the definite integral 

rX m ^{E) dx 

P = 2 - - = 42 

J-x m ME) V2^F{X) + E 

where X m3jX (E) is the positive root of F(X) + E = 0. The period P cannot be computed 
in closed form given the particular potential —F(X) in this problem. Alternately, one can 
compute the period simply by integrating equation (|4*T|) with initial conditions X = X max (E), 
X = until reaching the termination condition X = —X maiX (E) at the time P(E)/2. 

In either case, the above calculation must be inverted numerically to yield the energy as 
a function of the period, using the secant method or some variant. In the scaled variables, 
we have, from equation (fTUj) . 

P(E n ) = (2n + l)n 
2 A 

This is solved for E n = — and the resonant velocity is found using (jHSJ). This, and the 
other two approximations are shown in figure OH All solutions with positive n up to n = 14 
are shown in the figure. Approximations (|34|) and (}4T)|) both predict the existence of a 2-3 
window, while the numerical calculation ()43j) does not, and no such window is found by 
direct numerical simulation. 



E. Generalizaton to near-resonances and higher resonances 

The 2-pass resonant solutions are a countable, and thus measure-zero, family of initial 
conditions. Each 2-pass window has finite width whose left and right edges can be found by 
imposing the conditions that AE 2 = M 2 /2, so that the output energy is identically zero. It 
can be shown that the window widths scale as n~ 3 for large n. 

In between the 2-pass windows there is a complicated structure consisting of many nar- 
rower windows. These include 3-pass windows, which can be found as follows. A three-pass 
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FIG. 9: The resonant velocities, indexed by the number of complete oscillations of w(t), with 
(3 = 0.05. The thick solid curve at the bottom is the result of direct numerical simulation. From 
top to bottom, the other curves are the asymptotic results of equations and (|4Uj). and the 
value involving numerical calculation of the energy level, given the resonant period (|43[) 

resonant solution has three energy jumps. Just as W(t) and X(t) are even functions about 
t* in two-pass solutions, in 3-pass solutions, W(t) and X(t) are odd functions about their 
center time. We can place the three "center times" at t = —t , t = 0, and t = t , and notice 
that if the solution is odd, then AE = at t = 0. The change of energy at the second jump 
is AE = — 3/ c 2 00 /7T 2 (| + Xt ), which implies 



A corrected formula, as in equation (J47i|) . and a more accurate numerical condition, as in 
equation (|4Hjl . may also be derived. A general formula for the locations of higher complete 
resonances can be derived as in j^, but this equation must be solved numerically. 



(2n + l±|)7t 



A 



and gives 3-pass resonant solutions with 
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t 

FIG. 10: The solutions components of the 2-5 resonance of (a) 1 + W{t) of the simplified equa- 
tion © and (b) w(t) of the full ODE @. 

VI. THE EFFECT OF COUPLING TO A WEAKLY NONLINEAR OSCILLATOR 

We briefly discuss the discrepancies between the full ODE model and the simplified 
model ©, in order to account for the marked difference between the window locations 
between figures El and The most obvious, seen by comparing the two graphs of figure [TPl 
where the component w(t) and 1 + W(t) are both plotted for the 2-5 resonance with (3 = 0.2. 
It is clear that w(t), the fully nonlinear oscillation, has a larger amplitude, a larger period, 
and the mean about which it oscillates is displaced upward. As done following equation ()19j) . 
we fit t 2 — t\ from our numerical simulations and find t 2 — t\ w 5.00 (n + |) + 0.541 when 
(3 = 0.05. For the case (3 = 0.2, careful examination shows that t 2 — £i is not approximated 
that well by a linear fit, with the growth in t 2 — t\ slowing as n increased, and t 2 — t\ w 
6.38 (n + g) + 0.327 when the first 10 resonances are used, and t 2 — t\ ps 6.20 (n + |) + 1.93 
when resonances 11 through 20 are used, and the error in this fit is much larger than in the 
simplified model, especially when the leftmost windows are included. We see then that, in 
addition to a large correction to the oscillation frequency, a significant phase shift appears 
in the fully nonlinear dynamics. 

This shift in amplitude, frequency, and mean value can be explained using a strained 
coordinate or Poincare-Lindstedt expansion j^j . We assume that w = W + 1 satisfies the 
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homogeneous part of (jH 

d 2 w x 2 (\ 



1 (44) 



dt 2 w 2 \w 

where we have used scaling (jHJ) and have again set K = 1. Expanding this as a power series 
in W, we find 

d LK + X 2 w = X 2 (3W 2 - 6W 3 + 0{W A )) (45) 
We look for periodic solutions with initial conditions 

W(0) = e and W(0) = 0, (46) 

noting that these initial conditions are somewhat arbitrary, and leaving e positive and small, 
but for now undefined. We expand both the function W and the frequency of oscillations Q 
in powers of e 



W = J2e k W k . 1 (T) 

k=l 

oo 

n = J2 ekQ 



t 1 (47) 



where T = Qt, and assume that the solution has period 2tt in T. The equation is satisfied 
at each order in e, with the chosen to suppress secular growth terms. We find that 

n = A; O x = 0; fi 2 = -h 

and 

W = cos fit 

3 1 

Wi = - - cos nt - - cos 2£lt (48) 

13 3 

W 2 = -3 H cos fit + cos 2f2t + - cos 3fit 

8 8 

Thus the period of oscillation is decreased at larger amplitudes, as found from the least 
squares fits. 

It remains to determine a suitable value of e in (j4^j) and its effect on the resonance. The 
full ODE model (j3J) conserves the Hamiltonian 

2?r 2 . 2 , 8 , 1 2\ 16/3 „ /^X 



H=-X 2 + —w 2 + - 1 + — - - - — — F - , (49) 
2 3 3\ w 2 w w \ w ' 
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not written here in canonical variables. We see from figure |3] that v c , and hence AE is 
approximately the same for the full and simplified ODE systems. Expanding this in powers 
of W — w — 1, we obtain the approximate Hamiltonian 



if G(X)W is small enough to be ignored, by the solution to — \v\ = — 16/?F(X max ). 
Plugging this back into the energy (pf9*j) and using the expansion (jBHj) only for the coupling 
F(X/w) term, we obtain 



For resonant velocities sufficiently close to v c or for small enough /3, 16/5G(X max ) is negligibly 
small, and we can solve the resulting quadratic equation for W(0) and obtain 



For larger values of 16j3G(X mSuX ), the equation is cubic in W and the roots may be found by 
a perturbation expansion around the previously found roots. We use this value of W(0) as 
our value of e. For (3 = 0.05, it is sufficient to use e = y/3v c /4:, which gives a period 5.00, 
as was found from the linear fit. For f3 = 0.2, we find that weakly nonlinear theory is not 
useful as the first several terms of the expansion of the Hamiltonian in W of (|50|) are all 
found numerically to about the same order, so that ignoring the W 4 and W 5 term in the 
Poincare-Lindstedt expansion is invalid. 

VII. CONCLUSION 

We have explained many of the phenomena seen in the collision of vector solitons in CNLS. 
First we derived and justified a simplified version of a model derived by Ueda and Kath. 
Using a Melnikov integral, we estimated the critical velocity, and using matched asymptotic 
expansions near separatrices, we explained how to connect subsequent passes to construct 
an approximate solution. Imposing the condition that the total energy change after two 
passes is zero allowed us to find the locations of the exact two-pass resonant velocities, the 




As t — > — oo, all of the energy is stored as kinetic energy in the soliton modes H = ^Vq. At 
the symmetry time, which we can set to t* = 0, X = and w — 0, and X = X max is given, 




W(0) 
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centers of the two-pass windows. It remains to be seen if this phenomenon can be produced 
in physical experiments, but the experimental setup would appear to be fairly simple. 

More importantly, we have elucidated the mechanism underlying two-pass and two- 
bounce resonance phenomena in general. The important elements of an ODE model are 
the following: 

• a "position" mode Z(t) that moves in a potential well V(Z), which is localized near 
Z = 0, so that the force approaches zero at large distances, 

• a secondary oscillator mode W(t) that acts as a temporary energy reservoir; 

• a term C(Z,W) that couples the two modes together, also localized near Z = 0, so 
that coupling decays at large distances. 

Trapping takes place on the initial interaction if enough energy leaks from the position mode 
(Z) to the energy reservoir iW). In that case, the first mode crosses a separatrix curve in its 
unperturbed phase space. The energy change on subsequent interactions may be positive or 
negative, depending sensitively on the phase of W(t) at the interaction time, even though 
W may remain exponentially small. Eventually, in any Hamiltonian model, enough energy 
will eventually be transferred back to the position mode that it returns to the unbounded 
portions of phase space and escapes. It is transmitted if it escapes to +00 and reflected if it 
escapes to —00. All the models we have seen are Hamiltonian, but it should be possible to 
carry through much of the analysis in the presence of a simple dissipative term. 

We believe this mechanism to be present in all the systems which have displayed two- 
papers on kink-antikink interac- 
. A finite dimensional model for 
kink-antikink interactions in the 4 model is presented by Anninos et al. p. The model 
they derive is essentially of the form described, but the potential and coupling terms are 
much more complicated than (j^J), and it would be quite difficult even to find expansions 
about poles in the solution, as was done in section HVl of the current papers. What's more 
there is no natural small parameter measuring the coupling and the difference in time scales 
of the two modes, as we have here. The topology is slightly different in the 4 kink-antikink 
problem: the separatrix is given by an orbit homoclinic to infinity, rather than heteroclinic. 
This is indeed why the interaction is a "bounce" rather than a "pass." We have developed 



bounce resonance phenomena, including the foundationa 



tions in nonlinear Klein-Gordon equations 



y y y 
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a simple model that displays the same topology as in |5j and are in the process of analyzing 
it as a next step in understanding the two-bounce case. 
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